The programs/scripts (CentralityFlowThrough.py, KMC_SelectInitialSite.cpp, FindPaths.cpp) in this directory are adapted versions
of the programs used in:

K.S. Gomez-Haibach and M.A. Gomez: Revised Centrality Measures Tell a Robust Story of Ion Conduction in Solids.
The Journal of Physical Chemistry B 127, 9258 (2023).

The code for a single proton moving is used but there are commented out lines to give the reader an idea of what would need to 
be done if a site or vertex was removed as in the paper above.

**Need before running**
To work properly, the two C++ programs require nrutil.c and nrutil.h from Numerical Recipes in C second edition.This is available 
under "Obsolete Versions" in https://numerical.recipes/ - Look at Appendix B for both files if your team does not have a license.
To access these routines from this site, download the "References and Appendices" as copying and pasting the code directly from 
the book will affect the formatting. After downloading in PDF form, copy and pasting will not affect the formatting and the code 
can be moved to text files. The recipes are a great addition to any group. The kinetic Monte Carlo code also uses the random number 
generator ran3.cpp which is also available in the second edition of Numerical Recipes in C located in Chapter 7 "Random Numbers" 
on page 283. This code is given with comments which must be deleted in order for the code to run properly.The ran3 generator is 
still superior to the default random number generator for Monte Carlo and other numerical simulations.

To use a c code like nrutil.c within a c++ program, you need to declare the routines as external C in nrutil.h.  
https://en.cppreference.com/w/cpp/language/language_linkage.html shows you how to wrap declarations in the extern command.

**Data files used**

The three programs/scripts use data found in Energies.binding and Energies.TS

Energies.binding has the binding sites in the form 

Number x          y               z            Energy
1 0.50770549257 0.243852730271 0.115869991832 -331.91974

The x,y,z are in fraction of the box units and the energy is in eV

Energies.TS has the transition states in the form
from to type energy 
1 2 2 -331.644666 

The centrality script also uses the file POSCAR which is a VASP file with the positions of all the non-proton ions
in this example.

** How to run Flow Through Centrality**

The python3 script CentralityFlowThrough.py run with

python3 CentralityFlowThrough.py

reads the Energies.binding, Energies.TS, and POSCAR files and creates CentralityFlowThrough.vesta
which can be opened with the VESTA visualization program available at https://jp-minerals.org/vesta/en/ .

The script also creates a file Sites with the proton sites or site vertices in the form
number x y z Bolztmann Probability
a file called Connections 
fromSite toSite probabilityOfMove rateConstantOfMove
These use transition state theory with a prefactor of 1.0/ps to capture the time scale of the moves but not their
frequency differences so in effect are considering mainly the effect of energy in 
centrality.  Other prefactors could be read.

Note that this script and the other programs have variables like the system Temperature which you can change.

** How to run kMC **

The kmc program reads Energies.binding and Energies.TS and must be compiled in the following way
 gcc -c nrutil.c -o nrutil.o
 g++ KMC_SelectInitialSite.cpp nrutil.o -o kmc

To run the program, a kmc.in input file is provided which contains the lines

96  - number of binding sites
1000.0 - temperature in kelvin
-5871 - random number seed for the ran3.cpp script from numerical recipes
10000.0 - max time to run kmc for in ps or the unit determined by the rate constant prefactor in your data
33 - starting proton site

to run:
./kmc <kmc.in >kmc.out

The program is set up to do 1000 crossings (numberCrossings=1000;) but you can change it.  For each crossing, the highest 
barrier is recorded and the limiting barriers for each conducting trajectory are averaged. The output also allows you to
compare the fraction of the time spent at each site with the Boltzmann probability so you can check that the code is working.
To get an exact match, increase the number of crossings.  The kmc.out file prints out some of the input, verifies that the 
random number generator is working well, shows the crossing times and limiting barriers for individual trajectories, 
provides averages over trajectories, and lets you check that you are sampling correctly.  To visualize the region travelled,
the sites visited can be labelled in the vesta file found in the centrality program.

** How to run Find Paths**

The FindPaths program reads Energies.binding and Energies.TS and Centrality.out (from python2 centrality code above).
It is compiled in the following way
gcc -c nrutil.c -o nrutil.o
g++ FindPaths.cpp nrutil.o -o findpaths
but memory use could need to be updated depending on the size of path considered.

For the program, a find.in input file is provided which contains the lines

96  - number of binding sites
1000.0 - temperature in kelvin
1 - first binding site to consider
96 - last binding site to consider
8 - smallest number of steps in the path
8 - largest number of steps in the path
0.96 - fraction of the simulation box needed for a path to be considered as crossing the simulation box

to run:
./findpaths <find.in >find.out

The program outputs the average limiting barrier for all paths at the temperature considered as well as the fraction of 
these paths in the X, Y, and Z directions.  The max barriers in each direction are also recorded as well as the 
fraction of these limiting barriers that are rotations, intraoctedral transfers, and interoctahedral transfers.
The total number of paths found is given as well as the sites visited in the best path of the lengths considered.
If the number of paths is not too large, the user can comment out the exit below the print line:
"To see sorted paths, comment out exit in program below this print line in code", recompile and run again.
This will give details of every path found ordered from the most probable to the least probable. 
To visualize a path, the path sites can be labelled in the vesta file found in the centrality program.
